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Abstract 

The dynamics of a qubit under the decoherence of a two level fluctuator (TLF) in addition to 
its coupling to a bosonic bath is investigated theoretically. Two different methods are applied and 
compared for this problem. One is a perturbation method based on a unitary transformation. 
With the merit of our unitary transformation, non-adiabatic effect can be taken into account 
efficiently And the other one is the numerically exact method, namely the quasi-adiabatic path- 
integral (QUAPI) propagator technique. We find that the analytical method works well for a wide 
parameter range and show good agreement with QUAPI. On the other hand, The enhancement 
and the reduction of quantum decoherence of the qubit is checked with varying bath temperature 
T and TLF-bath coupling. 
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I. INTRODUCTION 



Dissipative quantum dynamics is one of the central paradigm in the theoretical physics 
because of its relation to the various physical and chemical phenomena from the spontaneous 
emission to electron transfer in molecular, from qubit decoherence to photon harvest in 

nn 

photosynthesis [1H5(. In the last decades, spin-boson model, the simplest possible model to 
describe dissipation, is studied intensively and offers a comprehensive understanding of the 
decoherence phenomenon. 

In this paper, a variant of the spin-boson model is studied, the overall spin-boson model 
act as an environment of another spin. This model describes a qubit in dissipative two- 



level fluctuators 
Josephson qubits 



F) which is believed to be relevant to the prevailing 1/f noise in the 



131 ] . Only a single spin-boson environment is considered since the qubit 
dynamics is usually dominated by a particular TLF near resonance with the qubit. An 
accurate evaluation of this problem is challenging. By employing the flow equation method, 
Gassmann et.al. studied the spectrum of the correlation functions of this model. And it is 
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-[l8|. 



studied a lot recently with various perturbation approaches 

In this paper we examine the effect of the bath temperature and system bath coupling on 
the qubit dynamics. One commonly believed concept is that temperature and the coupling 
to noisy bath only play a negative role in preserving the qubit coherence. However, it is 



pointed out in Ref. [19[ that the temperature can help the coherence when the qubit is 
coupled to a TLF (or spin-boson) environment. To examine this effect we treat the model 
more rigorously by consider spin-spin as a central quantum system which coupled to a 
bosonic bath and compared with the result of the numerically exact method, namely the 
quasi-adiabatic path-integral (QUAPI) propagator technique. Good agreement is achieve 
between these two methods. And we find that it is possible for the reduction of decoherence 
with increasing temperature or with increasing TL-bath coupling, which verifies the previous 
findings. 

The model is given as (h = 1): H = Ha + Hab + Hb with 

A a a n go a b 

H A = ~Y a z > H AB = a x > 

A a B 

H B+b ath = -Y^ + J2 Ukb l bk + ^Y,9k( b l + b k), (1) 

k k 

where TSSs are characterized by pseudospin-1/2 operators af and erf as usual, bk and b\ 



are the annihilation and creation operators of the bath mode, go and gk are the coupling 
constants. Here, the anisotropic coupling between TSS-A and TSS-B subject to the z di- 
rection, and the B-bath coupling also to the z direction. The bath is fully defined by the 
spectral density J{oj) = J2k9k^( u ~ u k)- We will use the piezoelectric spectral density, 
which describes the decoherence of a double quantum dots(DQD) qubit manufactured with 



GaAs 
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2l|, 



J 1L (u) =au[l 



— sm — 

u ui d 



(2) 



where Ud is related to the center to center distance, and ui to the dot size. Typically, 
Ud ~ O.Ol(ps)^ 1 and uj\ ~ l(ps) _1 (21|. In the limit of u d —> 0, one can find that, Eq. (J2J) 
goes back to the widely used Ohmic spectrum [22|, 123 ]. 



II. UNITARY TRANSFORMATION 



On the analogy with Ref. 24j, we apply a unitary transformation to the Hamiltonian, 
H' = exp(S)H exp(— S), with the generator S = J2k ~^k)o'x + i0o<Ty(T x 3 - Therefore, 



H' A 



(A^ cos 0q + g sin 9 ) 



h 'ab = x (9o cos e o - a a sin 0„) - ^A B sin 6 i(jyi<jy 



H 
H 
H 



1(0) 

B+bath 



B 



2 V A B cos 9 + "kblh + £^ (4 2 " 26 



B+bath 
,(2) 

B+bath 



T]A b COS 6*0 



AbC ° S9 ° \ a B (coshXi _ ) _ %0 B (sinhXi _ ^ 



+ 



A b sin 6q 



[erf sinhXi — i(Ty (coshXi — 77)] 



where X\ = J2k f^£fc(^l — ^k) an d V is the thermodynamic average of cosh A 1; 



rj = exp 



which insures H 



1(2) 

B+bath 



(3) 
(4) 
(5) 

(6) 



(7) 



(8) 



contains only the terms of two-boson and multi-boson non-diagonal 



transitions and its contribution to physical quantities is (gl) 2 and higher. Suppose the B- 
Bath coupling is not strong, the last term of the above Hamiltonian are dicarded in the 
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following discussion. Now if we let 

= cos 6*o + #o sin O , A' B = i]A B cosd , 
9o = (g cos9 - A A sm6 ) = r)A B sm9 , 



9k = 9k0- = r)A B cos9 Q — £ fc , 



then, 



tan O 



A A + r/A s 



and the Hamiltonian becomes 



H 



/(0) 



A' A' 

~Y z + ~Y z ' 



W fc + T^As COS O 



H 



/(0) 



r'(i) 

B+bath 



(9) 
(10) 

(11) 

(12) 
(13) 

(14) 
(15) 



where we have omitted the constant and the second order terms. 

It is easy to check that (h' ab + # B + bath ) \g ) = 0, where \g ) is the ground state of 
H' AB + Hfat h . Therefore, the ground state energy is E g ■ 



■ ±A' 



±A' 

2^B 



Note that, the transformed Hamiltonian H' 



r'(0) 
'bath 



Eft&(2-e*)- 



B+bath 



is of similar 



form of that of RWA, which enable us to treat the system bath coupling much easier, but 
A^, A B , g / 2 and g k /2 are replaced by A^, A' B and g' and gjj. due to the contributions 
of anti-rotating terms. The renormalized effective coupling become much smaller than the 
original coupling, which enable the pertubation treatment works better than the ordinary 
Born approximation directly from original Hamiltonian. 

Now, we can diagonalize H' A ^ B + H' AB by a unitary transformation T, 

10 

cos(0) sin(0) 

sin(0) -cos(0) 

1 

where tan 20 = 2g' /(A' A — A' B ). Then the Hamiltonian becomes H = H + Hi, 



T 



(16) 



H = ^2ei\i)(i\ + ^2uj k blb k 

i=0 k 

Hi = J29k(e+h + Q-bl)- 



(17) 
(18) 



where, e 3 = -e„ = ^ = \ (A' A + A' B ), e 2 = -e x = ^f 
cos#(|3>(2| - |1)(0|) + sin0(|2)(O| + |3)(1|). 



(A' A -A' B ) 2 + g'l e+ = e l 



III. MASTER EQUATION METHOD 



Now, we write out the master equation for this 4-level system by treating H as unper- 



turbed part, and Hi as perturbation, the reduced master equation is 

dp{t) 



25 



26|]: 



dt 



—i 



H e ,p(t) - / dt'X(t,t'). 



(19) 



where H e = X^=o £ *l*)(^l an< ^ X(t,t') is 



X(t, t) — TVbath 

2 n k e 



-iH () t 



Hup(t-t') 



JH Q t 



k 

+ 52 fin*-*** 



g_e- iH ^g + p(t-t')e lHEt - e~ iHet Q + p(t-t')e lHet q 



+ ^2g'k(n k + l)e- Wkt \g + e- iHet Q^p(t-t')e lHet - e~ iHet g^p(t-t')e lHst g+ 



iH,t 



The above master equation is a 4x4 matrix equation. According to the Kronecker product 
property and technique to Lyapunov matrix equation in matrix theory, by expanding the 
matrix p(t) into vector vec[p(t)] along row, the equation becomes, 



d 



H e ®l4 X 4 — hx4®H e 



^W)\ = ~J29k dt'F k (t-t')vec[p(t r )] 



with 



F k (t) = n k e™* % 
+ n k e~ i0J > 
+ (n k + l)e-™ kt 
+ (n k + l)e^ 



-iH e t 

f?- e Qa 



SH e t 



-iH e t 

e p_| 



T 



T 



e iH ^g+ 



T 



< (g+e^gY ' - ( g_e~ iS A ® ( g + e iS ^ '' 



(20) 
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In order to seek the non-Markovian effect of the non-linear bath, we solve the above 
equation by using Laplace transformation, rather than by replacing p(t') as p(t) which is the 
usual treatment in the literature known as the Markovian approximation. After the Laplace 
transformation and convolution theorem, the master equation of the system can be obtained 
as: 



U(P) 16xl6 vec p{P) 



vec 



(21) 



with 



U(P) 



16x16 



PL 



16x16 



+ i 



LL e ®L 4x4 - I 4x4 ®H e ] +Y^9k F k( p ) 



'16x16' 



and F k (P) 

i6xi6 ^ S ^ ne Laplace transformation of i*&(t)i6xi6- The master equation fl2T]) leads 
to following uncoupled equations: 



[P + iE p + (2 n k + 1) (sin 2 9B 7+ + sin 2 9B S+ )] p 14 (P) = pu(0) 
[P - iE p + (2n k + 1) (sin 2 9B 7 ^ + sin 2 0B 8 _)] p 41 (P) = p 4 i(0) 



.4 



44 " 



MP) MP) MP) MP) 
MP) MP) MP) MP) 
MP) MP) MP) MP) MP) MP) 



A' ■ 



Pi 2 (0) Pi 3 (0) p 24 (0) p 34 (0) 
P2l(0) p 3 i(0) p 42 (0) p 43 (0) 

1 T 



T 



(22) 
(23) 

(24) 

(25) 



1 T 



Pn(0) p 22 (0) p 23 (0) p 32 (0) p 33 (0) p 44 (0) 

(26) 

with the explicit definition of E> 7 ±, B$±, A 44 , A' M and A m being defined in Appendix. 

Before solving these equations, we first have a look at the initial condition and the physical 
quantities. In this work, the physcial quantity under our concern is the population difference 
Pit) = (<r£){t) = Tr [p(t) a £ (g>L 2x2] • It can be rewritten as, 



P(t) = Tr p(t)a£ 



where cx4 = Te s a A e S T = T(a A cos 9 + cr A a^ sin 9 )T, which is 

sm(9 + 9 ) -cos(9 + 9 ) 

sm(9 + 9 ) cos(9 + 9 ) 

-cos(9 + 9 ) sm(9 + 9 ) 

cos(9 + 9 ) sin(fl + fl ) 



(27) 



A 



(28) 
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thus, 



P(t) = {Pl2 + P34 + P21 + P43) Sin (0 + Bq) + (p 2 4 - Pl3 + P42 ~ P3l) COS (B + Bq) 

= 2Re(p l2 + p u ) sin (B + B ) + 2Re(p 24 - p 13 ) cos (B + B ) . (29) 

where the time dependant p(t) is replace by p for simplicity. From the above expression, we 
can see only Eq. f)24p is relevant to the dynamics of population difference. According to Eq. 
(124]) . we have 

(P + ds) [pi 2 (0) + p 34 (0)] + / [p 24 (0) - p 13 (0)] 



Pl 2 (P)+P 34 (P) 

p 24 (P) - Pi 3 (P) 



(P + rf 1 )(P + rf 3 )-/ 2 
(P + di) [p 24 (0) - p 13 (0)] + / [p 12 (0) + p 34 (0)] 



1+ 



(P + d 1 )(P + d 3 )-/ 2 

di = i(E p - E m ) /2 + ^2{2n k + 1) cos 2 BB 

k 

d 3 = i ( E p + P m ) /2 + ( 2n k + 1) sin 2 0Pi+ 
/ = ^ (2 n fc + 1) cos (B) sin (0) P 



'1+ 



(30) 
(31) 

(32) 
(33) 

(34) 



Suppose the system is in the upper eigenstate of erg and at the initial time t=0, therefore, 
the initial condition is given by 

r l 1 1 1 



P(0) 



1 


1 1 




1 1 


1 


1111 


2 


1 1 


<S>5 


1 1 


" 4 


1111 



1111 



(35) 



therefore, p(0) = |Te 5 (/ 2 4 x2 +^)(Jf x2 +af )e" 5 T = ±T(/£ 2 +<^ cos&o+afa* sin0 o )(/ 2 B x2 + 
&x)T, which is 

1 + sin (0 O ) cos (0) + sin (B + & ) sin (0) - cos (B + O ) cos (0 O ) 

cos (0) + sin (0 + 0o ) l + sin(20 + o ) -cos(20 + o ) cos (0 + O ) + sin (0) 

sin (0) - cos (0 + O ) -cos(2 + o ) l-sin(2 + o ) sin (0 + O ) - cos (0) 

cos (0 O ) cos (0 + O ) + sin (0) sin (0 + O ) - cos (0) 1 - sin (0 O ) 

from which we have 



P(0) = - 



Pi 2 (0)+p 34 (0) = sin(0 + o )/2 
p 24 (0)-pi 3 (0) = cos(0 + o )/2 



(36) 
(37) 



Therefore, the population difference is 

_! (P + d 3 ) sin 2 6 + / sin 6 cos 6 



P(t) = Re 



+ Re 



(P + d 1 ){P + d 3 )-P 
x (P + di) cos 2 9 + / sin 9 cos 9 



(P + d 1 )(P + 4)-/ 2 

where, J§f _1 is the Laplace inversion operator which corresponds to J„-i^ dPe pt , G 
9 + 9 and 



(38) 



o 



with 



(39) 
(40) 

(41) 
(42) 

(43) 

k uj + T]A B cos 9 J 

Here we would like to summarize the approximations we have made. Two approximations 
are made: The first one is the omission of H' 2 , which is a 4th order approximation to the B- 
bath coupling. The second one is the Born approximation for deriving the master equation 
( I19p . which is consistent with the first approximation. Therefore, our treatment is applicable 
in low temperature for a<Sl and g , T<^Aa, A s . 



di = i {E p - E m ) /2 + cos 2 9G{P) 
d 3 = i{E p + E m ) /2 + sin 2 9G(P) 



f = cos (9) sin (9) G(P) 

G(P)= / duo' ^iy cothOW), 
P + ico 

J {cu) = J[cu) 



IV. QUASI- ADIABATIC PATH-INTEGRAL METHOD 



The QUAPI method is a numerical scheme based on a exact methodology 27H30|. The 
starting point of QUAPI method is the generic system-bath Hamiltonian 

P 2 1 

1 * 2\2 



P- I 



2m 7 ' 2'- J " J ^ J 

3 J 

where, H is the Hamiltonian for the bare system, s is the system coordinate, and Qj 
are harmonic bath coordinates which are linearly coupled to the system coordinate. The 
characteristics of the bath are captured in the spectral density function 



The reduced density matrix of the system evolve as p(s", s', t) = 
Trbath (s"\e~ lHot p(0)e lHot \s'}. If the path integral representation is discretized by N time 
steps of length At = t/N and the initial density matrix is assumed to be p(0) = p s (Q)Pbath(ty, 
the reduced density matrix takes the form 

p( s », ^) = EE-EEEE e ~ lHoAt I 5 --) • • • W1 e ^ oAt 14 > 

S JV-1 S N-1 s l s l s S 

( s o \ Ps(0) I s ) (s I e |s 1 ) • • • (sjv-i| e lH ° At \s ) I(sq , s , sf , s x , • • • , s^_ 1 , s N _ 1 , s , s , At) 

(45) 

where the discrete variable representation (DVR) is used, the symbol (k = N — 1) 

denotes the system coordinate at the time kAt on the forward and backward discretized 
Feynman path. |s^) (k = iV — 1) are the eigenstates of the system coordinate op- 
erator s. If a symmetric splitting of the time-evolution operator is employed e~ lHAt = 
e -iH env At/2 e -iH At e iH env At/2 w ^ }j env — jj _ jj Qj the corresponding influence functional 

reads 

= TTh h \ e ~ iH env{s") At/2 e -iH env (s+_ 1 ) At e -iH env (s£) At/2 

X Pf, at h(0)e iHenv ( S ° } At l 2 ■ ■ ■ e iH env(s^_ 1 )At e iH env (s')At/2 ^ ^g-j 

One can find that the equilibrium position of the bath mode is adiabatically displaced 
along the system coordinate. If H provides a reasonable zeroth-order approximation to the 
dynamics, the quasi-adiabatic propagator is accurate for fairly large time steps. That is 
the quasi-adiabatic partitioning is a good representation when the bath property is mainly 
adiabatic, where the bath can keep up with the motion the system quickly. And the discrete 
path is to take into account of the non-adiabatic effect. For most of case, the quasi-adiabatic 
partitioning is reasonable especially when the system bath coupling is not strong. Therefore, 
the QUAPI discretization permits fairly large time steps when the adiabatic bath dominates 
the system dynamics. If the bath is purely adiabatic, even no discretization is needed. In the 
continuous limit (that is for At — > 0, N — > oo) the influence functional has been calculated 
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by Feynman and Vernon 

ft r f 



I = exp |-~ jf dt' dt" [s + {t') - s-(t')] 
x [a(t' - t")s + (t") - a*(t' - t")s-(t")] 



! l'*T4^ s+(f ' )2 - s " (f ' )2 ]| (47) 



where a(t) is the bath response function, which can be expressed in terms of the spectral 
density as 

i r°° r / R,,,..h \ i 

(48) 



a(t) = — I duJ(u) 



71 



(3tUjh 
2 



coth ( — ^— cos(cjjt) — ism(ojjt) 



The last term in Eg. (I47j) arises from the "counter-terms" which are grouped with the bath 
Hamiltonian in the quasi-adiabatic splitting of the propagator. With the quasi-adiabatic 
discretization of the path integral, the influence functional, Eq. H7J takes the form 



N k 




- s k] [vkk'4' -vlk' s k'] 

where = s" and = s'. The coefficients r] kk / can be obtained by substituting the 
discretized path into the Feynman- Vernon expression Eq. (]47l) . which is given in Ref. 28| . 

The QUAPI method is essentially a tensor multiplication scheme, which exploits the 
observation that for environments characterized by broad spectra the response function 
a(t) decays within a finite time interval. From the expression of the Feynman and Vernon 
influence funcitonal Eq. ( 145]) . one can see that a(t) characterizes nonlocal interactions, 
which connects system coordinate s(t') with s(t"). The path s ± (t / ) at time if is connected 
to the all the paths s~(t") at earlier times, which makes the evaluation of Eq. ( |45l) a hard 
task. However, for a bath with a broad spectral density, such as a power law distribution of 
the spectral density, a(t) has the finite memory, the memory length typically extending over 
only a few time slices when the quasi-adiabatic propagator is used to discretize the path 
integral. After discarding the negligible "long-distance interaction" with t' — t" > Ak max At 
(or k — k'> Ak max ), the resulting path integral can be evaluated iteratively by multiplication 
of a tensor of rank 2Ak max . In other words, there exists an augmented reduced density tensor 
of rank Ak max that obeys Markovian dynamics. The details of the multiplication scheme is 
discussed to a great extent in the literature, here we only present the essential parameters 



and mention briefly how to adopt it to our specific problem 27-3(3]. Here we discuss briefly 



the parameters used in the QUAPI method: 
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(i) The first parameter time-step At used for the quasi- adiabatic splitting of the path- 
integral. The memory time of the non-Markovian steps used by QUAPI is Ak max At. The 
stability of the iterative density matrix propagation ensures the choices of At, it should 
not be too big nor too small, since the non-adiabatic effect requires more splitting of the 
path integral, that is smaller At. Whereas, since the memory length Ak max At is usually 
a fixed value for a particular bath, QUAPI method prefers larger At, and consequently 
smaller Ak max in consideration of the numerical efficiency (note that the algorithm scales 
exponentially with A/c max , also see the discussion of the second parameter A/c max ). Therefore, 
we should choose appropriate At to take into account both the non-adiabatic effect which 
prefer smaller time splitting and the non-Markov effect which prefer long memory time, 
typically, we choose At around 2Q 2 ^ , that is to choose tens of fraction of the cycle time of 
the bare system dynamics. 

(ii) The second parameter is the memory steps Ak max . If Ak max < 1, the dynamics is 
purely Markovian. If the non-locality extends over longer time, terms with Ak max > 1 have 
to be included to obtain accurate results. In order to acquire converge result, in the practical 
implementation of QUAPI, one usually need to choose Afc max large enough so that the 
response function reduces to negligible value within the length of Ak max At. However this is 
a hard task, Since augmented propagator tensor A( Akmax > is a vector of dimension (M 2 ) Akmax 
(M is the system dimension which is four here), and the corresponding tensor propagator 
rp{2&k max ) j g a ma t r j x f dimension (M 2 ) 2Ahmax , the QUAPI scheme scales exponentially with 
the parameter A£; max . Thus one can not proceed the QUAPI calculation with very large 
AA; max , and usually A/c max is chosen less than 5 for M = 4, and even smaller for larger M. 

V. RESULTS AND DISCUSSION 

We report P(t) as a function of time in Fig. 1 to Fig. 4 for A^ = A# = O.lui, go = O.lA^ 
and Ud = 0.05a;;. The analytical results are depicted in solid lines, and QUAPI results of 
different Afc's are scatters with different colors and shapes. In Fig. 1, where a is larger than 
go/ A a which is set to be a = 0.3, it shows that the decoherence is reduced by increasing 
the bath temperature T as predicted in Ref. [19(]. However, in Fig. 2, where a = 0.01, the 
coherence is not meliorated but rather damaged with increasing T. Similarly, in Fig. 3, 
decoherence reduction with bath coupling when a = 0.3 is possible, whereas, in Fig. 4, there 
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is only decoherence enhancing where a = 0.01. 

One can see both analytical and numerical shows good agreement. When the coupling 
is small, e.g. a = 0.01,0.02 and 0.03 in Fig. 2 and Fig. 4, the two methods show perfect 
agreement, all the QUAPI results converges to our TRWA results as Ak increases. When a 
becomes larger, e.g. a = 0.3, 0.4 and 0.5, discrepancy appears, we attribute this discrepancy 
to that the QUAPI is not converged completely. Since when a becomes large the non- 
adiabatic boson contributes significantly, one need small At to take into account of the 
non-adiabatic bosons. Whereas, the memory length is almost the same, consequently, one 
can not converge within Ak = 3 which is the upper limit of our computation resources (note 
the algorithm scales exponentially as Ak). Finally, as a check of our analytical method, we 
report the strong qubit-TL coupling case, where go = A. One can see perfect agreement is 
achieved between the analytical and the numerical methods. 

VI. CONCLUSION 

In conclusion, without making RWA and Markov approximation, the dynamics of a qubit 
coupled with a spin-boson bath is investigated. We proposed an unitary transformation 
to treat this problem. The results of our analytical method show good agreement with 
QUAPI, even when the qubit-TL coupling is as strong as the TL spacing. And checked 
the decoherence behavior with varying bath temperature T and TLF-bath coupling. The 
decoherence of TSS-A can be reduced increasing bath temperature T or with increasing 
TL-bath coupling a only when the A-B coupling is smaller than B-bath coupling. 
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Appendix A: The matrix form and the matrix elements of ^fc(-P)i6xi6 anc ^ U(P)iq x 16 



From Eq. (1201) . the laplace transform of F^(t) is of following form, 
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(Al) 

where F k (P)ij are expressed as F i; j for simplicity. In order to save space, we do not give 



the explicit expressions of Fij here. On the other hand, H e ®I 
diagonal matrix. Therefore, from the definition of 6 r (P)i6xi6 ) 



4x4 



'4x4^ 



)H P 



is just a 



U(P) 



16x16 



-P-?16xl6 + i 



4x4 



16x16' 



(A2) 



we can get each element of < 7 7(P)i6xi6 as, 



U, 



m,n 



E k 9'k 2 F k (P)r 



when m ^ n 



Ur, 



P + i 



4x4 J 4x4 l < 



Y. k g' k 2 F k {P) 1 II, T 



With the particular form of U(P)iq x iq one can decouple the master equation Eq. (12T1) into 
equation sets with smaller dimension as shown in Eq. (}22l) - (l26l) . which are explicitly ex- 
pressed in the subsequenct Appendix. The parameters B±'s which will be used in the 
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master equation sets are defined as, 

1 



B 1± 
B 7 ± 



P ± iu k 



1 



2± 



1 



p±i(u k -E m y 
i 



p ± *k + 
i 

P± l (cu k -^fi) 



B, 



6± 



P±*(w fc + ^ 



Er, — E m \ ' 



5; 



8± 



1 



P±t(co k +^±^) 



Appendix B: Solve the 4x4 Master equation 

The master equation for pi 2 , pi 3 , p 24 and p3 4 is 

T 



Pu(P) Pis(P) P2a(P) MP)] = [pi 2 (0) pis(O) p 24 (0) p 34 (0)] 



where A 44 is defined as 



E,^ 1 ^ P + 4 + n^ 4fc 



2 



-n k e lk 
n k di k 



—n k d,2 k 
n k e lk 



(i) 

-n\ 'e lk 

— % Gt 2 £; 



(i) 
n k e lk 



-Ek^T^ P + d 1 + n k d 



2k 



where double indexes k indicate a sum over k, = n k + 1 and 

d 1 =i(E p - E m ) /2 + J2 9k (2 n k + 1) cos 2 9B 1+ 

k 

d 3 = i(E p +E m )/2 + J2 9k (2 n k + 1) sin 2 6B 1+ 



d 2k = g' k 2 sm 2 9(B 2+ + B 4 _) 

d ik = g' k 2 cos 2 6 (£3+ + #4-) 

ei* = ^ 2 cos(0)sin(0)(P 4 _ + P 1+ ) 

e 2fc = ^ 2 (2 n k + 1) cos (0) sin (9) (P 4 _ - P 1+ ) 

By solving the above 4x4 linear algebra equation, we can get 

(P + d 3 ) [p 12 (0) + p 34 (0)] + / [p 24 (0) - p 13 (0)] 



Pl 2 (P)+P34(P) 

MP) -MP) 



(p + rf 1 )(p + 4)-/ 2 

(P + jl) [P24(0) - Pl 3 (0)] + / [Pl2(0) + P34(0)] 
(P + rf 1 )(P + rf 3 )-/ 2 
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/ = [( 2 n k + 1) eik ~ e 2k ] /2 = ^ ^ 2 (2 n fe + 1) cos (0) sin (9) B 1+ 

k k 

As for the other 4x4 master equation, since p 2 i, P31, P42 and p 43 is just the complex 
conjugate of pi 2 , pi 3 , p 24 and p 34 , we do not need to solve them, and ,on the other hand, 
A' 44 is the very similar to A 44 , the only difference is that 5 + 's are changed to BJs and vice 
versa. 



Appendix C: Solve the 6x6 Master equation 



A 



The master equation for p n , p 22 , P23, P32, P33 and p 44 is 

1 T 



66' 



MP) MP) MP) MP) MP) MP) 



Pn(0) p 22 (0) p 23 (0) p 32 (0) p 33 (0) p 44 (0) 

(Cl) 



A 66 defined as 



P + n^ 1} (aifc + a 2k ) 


~ n k a lk 




-n k b 2 k 


—nkd 2 k 





(i) 


P + ^fcaifc + nj^W 


61+63 
2 


62+64 
2 





— nka 2 k 




2 


P + a 3 





61—63 
2 


nkbik 


-n k 1) b 2k 


62+64 
2 





P + a 4 


62 — 64 
2 


n k b 2k 


(i) 
; a 2fe 





61—63 
2 


62— 64 
2 


P + nf 1 l a lk + n k a 2k 


—nkCiik 





(i) 

-n\ 'a 2k 


n k ] b lk 


n k ] b 2k 


(i) 

-n k 'a^ 


P + n k (aifc + a 2k 
(C2) 



a 2 fc 
a 3 



= ?j 



where 61 = X fe &ifc, & 2 = Efe & 2fc, 63 = Efc & 3fc, &4 = Efc fe 4fc, 

^cos 2 0(S 6+ + S 6 _) 

^sill 2 0(S5 + + S5_) 

+ ^ 2 (2n fe + 1) (cos 2 9B^ + sin 2 #5 6+ ) 

k 

a 4 = -iE m + 9k ( 2n k + 1) (sin 2 9B^ + cos 2 #£5+) 

k 

b lk = g' k 2 cos(9)sm(9)(B 5 _ + B 6+ ) 
b 2k = g' k 2 cos(9)sm(9)(B 6 _ + Br 0+ ) 
hk = 9k (2n fc + 1) cos (9) sin (9) (5 5 _ - B 6+ ) 
Wk = 9k (2n fc + 1) cos (9) sin (6) (5 6 _ - B 5+ ) 



(C3) 
(C4) 
(C5) 

(C6) 

(C7) 
(C8) 
(C9) 
(CIO) 
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The analytical solution of this equation set is still out of our capability, however, it is 
irrelevant to the quantities which we are interested in. If we add the first two and last two 
lines of the master equation, we can find 

MP) + MP) + MP) + MP) = VP (en) 

which ensures the conservation of the trace of the reduced density operator Trp = 1. 



[1] Y. Nakamura, Y. A. Pashkin, and J. S. Tsai, Nature 398, 786 (1999). 

[2] I. Chiorescu, P. Bertet, K. Semba, Y. Nakamura, C. Harmans, and J. E. Mooij, Nature 431, 



159 (2004), URL http://dx.doi.org/10.1038/nature02831. 



[3] J. Q. You and F. Nori, Phys. Today 58, 42 (2005). 

[4] G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ami, T. Mancal, Y. C. 
Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007), URL 



http : //dx . doi . org/10 . 1038/nature05678, 



[5] H. Lee, Y. C. Cheng, and G. R. Fleming, Science 316, 1462 (2007), URL 

|http://dx. doi. org/10. 1126/science. 1142188] 
[6] E. Paladino, L. Faoro, G. Falci, and R. Fazio, Phys. Rev. Lett. 88, 228304 (2002), URL 

|http : //dx . doi . org/ 10 . 1 103/PhysRevLett . 88 . 228 304, 



[7] R. W. Simmonds, K. M. Lang, D. A. Hite, S. Nam, D. P. Pappas, and J. M. Martinis, Phys. 



Rev. Lett. 93, 077003 (2004), URL http://dx.doi.org/10.1103/PhysRevLett.93.077003 



[8] G. Falci, A. D'Arrigo, A. Mastellone, and E. Paladino, Phys. Rev. Lett. 94, 167002 (2005), 



URL http : //dx . doi . org/10 . 1 103/PhysRevLett . 94 . 167002 



[9] A. Shnirman, G. Schon, I. Martin, and Y. Makhlin, Phys. Rev. Lett. 94, 127002 (2005), URL 



htt p : //dx . doi . org/ 10 . 1 1 03/PhysRevLett .94. 12 7002, 
[10] B. Abel and F. Marquardt, Phys. Rev. B 78, 201302 (2008), URL 

|http : //dx . doi . org/ 10 . 1 103/PhysRevB . 78 . 201302| 
[11] M. Neeley, M. Ansmann, R. C. Bialczak, M. Hofheinz, N. Katz, E. Lucero, A. O'Connell, 

H. Wang, A. N. Cleland, and J. M. Martinis, Nat. Phys. 4, 523 (2008), URL 

|http : //dx . doi . org/10 . 1038/nphys972[ 



16 



[12] A. Lupa§cu, P. Bertet, E. F. C. Driessen, C. Harmans, and J. E. Mooij, Phys. Rev. B 80, 

172506 (2009), URL |http : //dx . doi . org/10 . 1103/PhysRevB . 80 . 172506| 
[13] J. Lisenfeld, C. Muller, J. H. Cole, P. Bushev, A. Lukashenko, A. Shnirman, and A. V. Ustinov, 



Phys. Rev. B 81, 100511 (2010), URL http://dx.doi.org/10.1103/PhysRevB.81.100511 



[14] H. Gassmann, F. Marquardt, and C. Bruder, Phys. Rev. E 66, 041111 (2002), URL 
|http : //dx . doi . org/10 . 1 103/PhysRevE . 66 . 041 1 11 



[15] H. Gassmann and C. Bruder, Phys. Rev. B 72, 035102 (2005), URL 
|http : //dx . doi . org/10 . 1103/PhysRevB . 72 . 035102 



[16] E. Paladino, M. Sassetti, G. Falci, and U. Weiss, Chem. Phys. 322, 98 (2006), URL 



|http : //dx . doi . org/10 . 1016/ j . chemphys . 2005 . 08 . 065[ 
[17] E. Paladino, M. Sassetti, G. Falci, and U. Weiss, Phys. Rev. B 77, 041303 (2008), URL 

|http : //dx . doi . org/10 . 1103/PhysRevB . 77 . 041303| 
[18] N. P. Oxtoby, A. Rivas, S. F. Huelga, and R. Fazio, New J. Phys. 11, 063028 (2009), URL 

|http: //dx. doi . org/10 . 1088/1367-2630/11/6/063028 



[19] A. Montina and F. T. Arecchi, Phys. Rev. Lett. 100, 120401 (2008), URL 



http : //dx . doi . org/10 . 110 3/PhysRevL ett . 100 . 12 0401 



[20] T. Brandes and B. Kramer, Phys. Rev. Lett. 83, 3021 (1999), URL 
|http : //dx . doi . org/10 . 1103/PhysRevLett . 83 . 3021 



[21] Z. J. Wu, K. D. Zhu, X. Z. Yuan, Y. W. Jiang, and H. Zheng, Phys. Rev. B 71, 205323 (2005), 



URL http : //dx . doi . org/10 . 1103/PhysRevB . 71 . 205323, 



[22] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. 



Mod. Phys. 59, 1 (1987), URL http://dx.doi.org/10.1103/RevModPhys.59.! 



[23] U. Weiss, Quantum Dissipative Systems (World Scientific, Singapore, 1999), 2nd ed. 

[24] H. Zheng, Eur. Phys. J. B 38, 559 (2004), URL 

|http : //dx . doi . org/10 . 1140/ep jb/e2004-00152-"7j 
[25] D. P. DiVincenzo and D. Loss, Phys. Rev. B 71, 035318 (2005), URL 

|http : //dx . doi . org/10 . 1103/PhysRevB . 71 . 035318 



[26] G. Burkard, Phys. Rev. B 79, 125317 (2009), URL 

|http : //dx . doi . org/10 . 1103/PhysRevB . 79 . 125317 



[27] D. E. Makarov and N. Makri, Chem. Phys. Lett. 221, 482 (1994), URL 
|http : //dx . doi . org/10 . 1016/0009-2614 (94) 00275^4 



17 



[28] N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4600 (1995), URL 

|http : //dx . doi . org/10 . 1063/1 . 469508] 
[29] N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4611 (1995), URL 

|http : //dx . doi . org/10 . 1063/1 . 469509 
[30] M. Thorwart, E. Paladino, and M. Grifoni, Chem. Phys. 296, 333 (2004), URL 

|http : //dx . doi . org/10 . 1016/ j . chemphys . 2003 . 10 . 007[ 

Figures Captions 

Fig. 1: P(t) is plotted as a function of time for different temperatures T/Aa = 1,5, 10 
with analytical method (solid lines) and the QUAPI methods with Ak = 1,2 and 3 
(respectively black square, red circle and green triangle). The decoherence is reduced 
with temperature T. The parameters: A^ = A_b, go — 0.2A^, a = 0.3, and the QUAPI 
parameter At = 0.4/A^. 

Fig. 2: P(t) is plotted as a function of time different temperatures T / Aa = 0.1, 1, 5 with 
analytical method (solid lines) and the QUAPI methods with Ak = 1,2 and 3 (respectively 
black square, red circle and green triangle). The decoherence is enhanced with T. Two 
frequencies are dominating the dynamics which agree with the results of Jaynes-Cummings 
model. The parameters: A^ = A B , g = O^A^, a = 0.01, and the QUAPI parameter 
At = 0.4/Aa. 

Fig. 3: P(t) is plotted as a function of time for different TL-bath coupling a = 0.3, 0.4 
and 0.5 with analytical method (solid lines) and the QUAPI methods with Ak = 1,2 and 3 
(respectively black square, red circle and green triangle). The decoherence is reduced with 
temperature T. The parameters: A^ = A B , g = 0.2A^, T = O.lA^, and the QUAPI 
parameter At = 0.4/ A^. 

Fig. 4: P(t) is plotted as a function of time for different TL-bath coupling a = 0.01, 0.02 
and 0.03 with analytical method (solid lines) and the QUAPI methods with Ak = 1,2 and 
3 (respectively black square, red circle and green triangle). The decoherence is reduced with 
temperature T. The parameters: A^ = A#, go = 0.2A^, T = O.IAa, and the QUAPI 
parameter At = 0.4/A^. 
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Fig. 5: P(t) is plotted as a function of time for strong qubit-TL coupling case g = 
Aa with analytical method (solid lines) and the QUAPI methods with Ak = 1,2 and 3 
(respectively black square, red circle and green triangle). The parameters: A^ = A B , 
a = 0.01, and the QUAPI parameter At = 0.4/A A . 
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